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Recently, realistic lattice QCD calculations with 2+1 flavors of domain wall fermions and the 
Iwasaki gauge action have been performed by the RBC and UKQCD collaborations. Here, results 
for the bottomonium spectrum computed on their gauge configurations of size 24 3 x 64 with a lattice 
spacing of approximately 0.11 fm and four different values for the light quark mass are presented. 
Improved lattice NRQCD is used to treat the b quarks inside the bottomonium. The results for the 
radial and orbital energy splittings are found to be in good agreement with experimental measure- 
ments, indicating that systematic errors are small. The calculation of the T(2S f ) — T(1S) energy 
splitting provides an independent determination of the lattice spacing. For the most physical en- 
semble it is found to be a -1 = 1.740(25)(19) GeV, where the first error is statistical/fitting and the 
second error is an estimate of the systematic errors due to the lattice NRQCD action. 

PACS numbers: 12.38.Gc, 14.40.Gx 



I. INTRODUCTION 

Bottomonium mesons, the bound states of bottom 
quark-antiquark pairs, play an important role in the 
study of the strong interactions. The spectrum of bot- 
tomonium is known very well from experiment and there 
are many different approaches to calculating it theoret- 
ically. Lattice QCD provides a model-independent and 
accurate way of doing this. 

One of the most important steps toward realistic lat- 
tice QCD calculations was the inclusion of dynamical 
light (u, d and s) sea quarks. For many quantities of 
phenomenological interest, lattice QCD now allows reli- 
able non-perturbative calculations that were previously 
impossible. In order to further control systematic errors 
and increase the confidence in lattice results, it is crucial 
to consider several different lattice actions and thereby 
test universality. 

The RBC and UKQCD collaborations have recently 
started large-scale lattice QCD calculations [l|, 0, H, 
EL S @ with dynamical domain wall fermions and 
renormalization-group- improved gauge actions. The do- 
main wall fermion action 0, d, @] has an approximate 
chiral symmetry that becomes exact, even at finite lat- 
tice spacing, when the extent L s of the auxiliary fifth 
dimension is taken to infinity. This leads to better con- 
trol over the renormalization of operators, a reduction of 
discretization errors and more reliable chiral extrapola- 
tions. 

The gauge configurations created by the RBC and 
UKQCD collaborations have been made publicly avail- 
able. In this work, the ensembles of size 24 3 x 64, L s = 16 
described in detail in [J] were used to compute the spec- 
trum of bottomonium. Chiral symmetry is not as im- 



portant for bottomonium as it is for light hadrons, but 
the calculations presented here nevertheless provide use- 
ful tests of the recent lattice calculations by the RBC 
and UKQCD collaborations. In particular, they provide 
an independent determination of the lattice spacing and 
a good way of tuning the b quark mass. The value of 
the bare b quark mass obtained here can also be used in 
lattice QCD calculations for heavy-light hadrons such as 
B mesons, which are of considerable importance for the 
phenomenology of weak decays and tests of the Standard 
Model. 

The b quark has a mass larger than the inverse lattice 
spacing, so that the standard lattice actions as used for 
light quarks are not suitable to describe it. A prelim- 
inary calculation of the bottomonium spectrum on the 
RBC/UKQCD gauge configurations using a relativistic 
heavy-quark action was presented in llCW . Here, improved 
non-relativistic QCD (NRQCD) [Hill [3 is employed 
instead, following closely the methods used in the calcu- 
lation of the bottomonium spectrum on the MILC gauge 
configurations in [14| (these configurations use an im- 
proved staggered fermion action for the sea quarks and a 
one-loop Symanzik- improved gluon action). 

The lattice calculation and analysis methods are de- 
scribed in detail in Sec. |TTJ The results for the tuning of 
the b quark mass and tests of the dispersion relation are 
given in Sec. IIII Al and IIIIB1 respectively. In Sec. IIII CI 
the results for the radial and orbital energy splittings as 
well as the lattice spacing are presented, followed by the 
fine and hyperfine structure in Sec. IIII Dt 



II. THE LATTICE CALCULATION 



A. Lattice actions and parameters 
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The details of the domain wall fermion and Iwasaki 
gauge actions as used by the RBC and UKQCD collab- 
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ami 


UOL 


MD range (step) 


^conf 


0.005 


0.8439 


915 


- 8665 (25) 


311 


0.01 


0.8439 


1475 


- 8525 (25) 


283 


0.02 


0.8433 


1800 


- 3600 (25) 


73 


0.03 


0.8428 


1275 


- 3050 (25) 


72 



TABLE I: The ensembles of RBC/UKQCD gauge configura- 
tions used here. 

orations are given in [lj]. Here, thegauge configurations 
of size 24 3 x 64 as described in [J] were used. These 
have L s = 16, f) = 2.13 and the strange quark mass 
is am s = 0.04. There are ensembles with four differ- 
ent values for the degenerate light (up and down) quark 
mass am;, as shown in Table |TJ The "measurements" in 
this work were started at the same molecular dynamics 
(MD) time as in 0] to ensure complete thermalization. 
Note however that the ami = 0.005 and ami — 0.01 
ensembles have since been extended and the additional 
configurations were included here. Measurements were 
performed every 25 steps of MD time, as discussed fur- 
ther in Sec. PTUl 

The lattice NRQCD action for the b quark is the same 
as in the previous study of the bottomonium spectrum 
in [ijj], with stability parameter n = 2. The full details 
of the action can be found in e.g. [l5j]. After the initial 
tuning, which will be described in Sec. IIII1 the bare b 
quark mass was set to ami, = 2.536. 

All couplings in the NRQCD action are set to their 
tree-level values, but the action is tadpole-improved 
which accounts for a large amount of the renormaliza- 
tion. The mean link in Landau gauge uql is used as the 
tadpole improvement parameter; the values of uol for 
the different ensembles are listed in Table HI 

The NRQCD action in use includes relativistic cor- 
rection terms up to order 0(v ) where v is the internal 
speed of the b quarks inside the meson. For bottomo- 
nium, one has v 4 ~ 0.01. The action is also tree-level 
Symanzik improved. Systematic errors depend strongly 
on the observable under consideration, and will be dis- 
cussed individually in Sec. IIIII However, finite-volume 
errors are expected to be negligible in all cases due to 
the small size of the bottomonium mesons. 



B. Calculation of the meson two-point functions 

As in [l3[ , the heavy- heavy meson correlators with mo- 
mentum p were computed from 

c(T sk ,r sc , P ,t-t r ) 

= Tr [&(x!,t, x[ 1 t')Tl ]i (x 1 -x 2 )G sc (x 2l t 1 x[,t') 

331,052 

x e-^-^ (1) 



where t > t' and 

G sc {x 2l t, x'^t') = ^G(:r 2 ,<, x' 2 ,t')T sc (x[- x' 2 ) 

xe lp ^-^. (2) 

In Eqs. Jl} and which are understood to be for a 
single gauge configuration, G denotes the heavy-quark 
propagator which is 2 x 2 matrix-valued in spinor space 
and 3x3 matrix- valued in color space. The functions 
r sc /sk are the "smearing functions" at source and sink, 
respectively, which are also 2x2 matrix-valued in spinor 
space. No gauge links were included in r sc / sk ; instead 
the gauge configurations were fixed to Coulomb gauge. 

In Table [Til the bottomonium states considered in this 
work are listed, together with their continuum quantum 
numbers, smearing functions T(r) and representations of 
the octahedral group As can be seen in the table, 

all representations are chosen to be different, so that no 
mixing is expected here. The radial functions </) n s(T), 
4>np(x) and 4>nD{f) for the n-th radially excited 5- wave 
[L = 0), P-wave (L = 1) and D-wave (L = 2) states 
were taken from the corresponding hydrogen atom wave 
functions and are given in Table IIIII The same lattice 
representations were used at source and sink but the ra- 
dial smearing functions were allowed to be different. The 
smearing parameters ro (in lattice units) were set to 1.0 
(15), 0.8 (25), 0.6 (35), 0.5 (IP), 0.4 (2P) and 0.5 (ID), 
respectively. 

Note that G in Eq. © can be computed efficiently by 
using the function 

r^aci-x^e*^ (3) 

as the initial condition in the heavy quark evolution equa- 
tion. For some states it is computationally more conve- 
nient to remove the Pauli matrix in r sc in the initial 
condition, and instead including it explicitly in the trace 
in Eq. (jl}. In this way, different spin directions can be 
obtained with a single G. 

Finally, note that on a finite lattice, the smearing func- 
tions must satisfy the periodic boundary conditions. This 
was ensured by setting the smearing functions to zero 
outside a ball with radius R smaller than half the spatial 
lattice dimension. In this way, the wrapping around the 
lattice boundaries does not cause any problems. Since 
the smearing functions decay exponentially with the sep- 
aration between quark and antiquark, R can be chosen 
such that the important features remain. To ensure sym- 
metry, the same cut-off radius must be taken at source 
and sink. 

In order to increase statistics, the correlators ([1]) were 
averaged over eight different spatial origins x[ located at 
the corners of a cube with side length L/2 = 12. In ad- 
dition, four different source time slices t' with an equal 
spacing of 16 were used, thus leading to a total of 32 ori- 
gins per configuration. Furthermore, the locations of the 
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TABLE II: The smearing functions T(r). See e.g. [l7|| for the irreducible representations of the octahedral group. 



State 



IS exp[-|r|/ro] 

25 [1 - |r[/(2ro)] exp[-|r|/(2r )] 

3S [1 -2|r|/(3r ) + 2|r] 2 /(27rg)] exp[-[r|/(3r )] 

IP exp[-|r|/(2r )] 

2P [1- |r|/(6ro)]exp[-|r|/(3ro)] 

ID exp[-|r|/(3r )] 

TABLE III: The radial functions 6(r) 



origins on the lattice were shifted randomly from config- 
uration to configuration in order to decrease autocorre- 
lations. 



C. Fitting and analysis details 

1. Bayesian multi- exponential fitting 

After choosing a set of smearing functions T(r) with 
equal lattice representations but different radial functions 
(e.g. IS, 2S and 35), the square matrix of correlators 
obtained by taking all combinations for source and sink 
was computed. 

The matrix of correlators (C(T s k, r sc ,p, t — t')) was 
simultaneously fitted by a function of the form 



n=0 



A n (T sc ) A* n (T sk ) e -^(*-*') 



(4) 



where E n is the energy of n-th state and A n (T) is the 
(real) amplitude for this state to be created by the oper- 
ator with smearing function T(r). 

To ensure the correct ordering of the states in terms 
of their energy, the fit parameters were actually chosen 
to be the logarithms of the energy differences between 
neighboring states (in lattice units) 



and the logarithm of the ground state energy, ln(i?o). 
Furthermore, the amplitudes for the excited states were 
written as 



A n (T) = A' n (T)A (T), 



(6) 



taking the relative amplitudes A' n (T) (for n > 1) and the 
ground state amplitude ^4o(r) as the fit parameters. 

The Bayesian fitting method described in [l8| was 
used, where the \ 2 function is augmented by 



X 



with the Gaussian prior 



X +X 



x 2 ■ 

A prior 



E 



ipi 



prior 



(7) 



(8) 



HE, 



n+l 



E n ) 



(5) 



Here, { Pi } = {A {T), A' n (T), ln(£ ), \n(E n+1 -E n )} are 
the fitting parameters, and the prior for each parameter 
Pi is given by its central value pi and width ap i . 

The Bayesian method allows the inclusion of an arbi- 
trary number of exponentials n eX p in © and hence the 
fitting in the full range of Euclidean time t — t' between 
source and sink. Here, only the points with t — t' = 
were excluded in the fits. The number of exponentials is 
increased until the fit results and error estimates become 
independent of n eX p- This is demonstrated in Fig. [T] 

In the following discussion of the prior choices {pi , ap i } 
we will distinguish between parameters for low-lying and 
high-lying states. For example, in a 3 x 3 matrix fit 
containing sources optimized for the T(15), T(25) and 
T(35) states, these three states will be referred to as 
low-lying, as their energies and amplitudes will be well- 
determined by the data, while higher excitations will be 
referred to as high-lying. 

The prior widths for the parameters of the low-lying 
states were chosen to be about 10 times larger than the 
resulting error estimates from the fit. This ensures that 
the influence of the priors on these parameters is negli- 
gible. Initial guesses for the central values were obtained 
from unconstrained fits including only a small number of 
exponentials at large Euclidean time. 
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FIG. 1: Fit results for the 3x3 matrix correlator with the 
{T(IS'), T(2S), T(3S')} smearings as a function of the num- 
ber of exponentials. The values of x 2 divided by the number 
of degrees of freedom are also shown. The results are for the 
ensemble with ami = 0.005. 



For the high-lying states, the priors for the loga- 
rithms of the energy-splittings between successive states, 
ln(E n+ i — E n ), were set to —1.4 (corresponding to about 
400 MeV) with a width of 1. The priors for the relative 
amplitudes A' n (T) of the high- lying states were set to zero 
with a width of 5. 



2. NRQCD energy shift 

Due to the use of lattice NRQCD, where the heavy 
quark mass has been integrated out of the action, all 
energies obtained from the fits are shifted, 



E : 



v phys 



2C 



(9) 



where C is approximately equal to the heavy quark mass. 
Since C is the same for all states, energy splittings are 
unaffected by this shift. 

The physical mass of a meson (and hence C) can be 
calculated from the energy difference between its states 
with p = and p ^ 0. Assuming the relativistic contin- 
uum dispersion relation 



E(p) = vV + M 2 - 2 C 
we obtain the kinetic mass 

p 2 - [E(p) - E(Q)} 2 



M = M v 



2 [E(p) - £7(0)] 



(10) 



(11) 



With the improved lattice actions used in this work, the 
continuum dispersion relation (|10p was found to be an 
excellent approximation for bottomonium at small lattice 
momenta p. This will be demonstrated in Sec. IIII Bl 



3. Bootstrap method 

When computing quantities that depend on more than 
one fit result, such as the kinetic mass (ITU or an energy 
splitting obtained from independent fits, correlations be- 
tween the different fit results must be taken into account. 
In this work, the bootstrap method was employed to 
achieve this. Note however that it must be modified for 
Bayesian fitting [l8[ so that not only the data sets are 
resampled randomly but also the central values pi of the 
priors in ([8]) are drawn from Gaussian random distribu- 
tions with widths ap i for every fit. 

The final quantity of interest is then computed for the 
bootstrap ensemble of fit results, giving an approximate 
probability distribution. In the end, the mean value and 
the 68% width of this distribution are quoted. 

The bootstrap method was in fact not only used for 
quantities depending on more than one fit parameter but 
also to obtain the error estimates for individual fit pa- 
rameters. The number of bootstrap samples was taken 
to be 500. 



4- Autocorrelations 

The integrated autocorrelation time Tj nt for the 12th 
time slice of the pion correlator on the ami = 0.005 en- 
semble was found to be 10 to 15 steps of molecular dy- 
namics (MD) time in (4j. Therefore only gauge config- 
urations separated by 25 steps, which is approximately 
equal to 2Tj nt , were used here. However, the autocorrela- 
tion time depends on the observable. The measurements 
in this work were checked for residual autocorrelations 
using the binning method. 

Recall from Sec. Ill Bl that on each gauge configuration 
meson correlators from 32 origins (8 origins on 4 source 
time slices each) were computed. The data was always 
averaged over the 8 origins on each source time slice, thus 
leaving 4 data samples per configuration. 

To estimate the autocorrelations between the different 
gauge configurations, the data was also averaged over 
the 4 source time slices prior to the binning. Note that 
the measurements already had an initial separation of 25 
steps in MD time, and hence the binning increases this 
to integer multiples of 25. 

Of course the binning reduces the number of data sam- 
ples available for the fit. To obtain a reliable estimate of 
the data covariance matrix (see e.g. [181]), the number of 
data samples should be much larger than the dimension 
of this matrix. Thus, in the analysis of autocorrelations, 
fits with only one smearing at source and sink and a small 
fitting range, corresponding to a small data covariance 
matrix, were considered. 

No significant increases in the bootstrap errors were 
seen for any of the ensembles, indicating that the sepa- 
ration of 25 steps of MD time gives sufficiently indepen- 
dent measurements. Note that the origins of the meson 
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correlators were shifted randomly between gauge config- 
urations. 

Next, tests for autocorrelations between the data sam- 
ples from the four different source time slices were per- 
formed. To this end the data was averaged into bins of 
2 or 4 time slices, without additional binning over gauge 
configurations. Here, in some cases a slight increase in 
the bootstrap errors was seen, at most 20%. Thus, for 
the measurements in the remainder of this work the fol- 
lowing conventions were used: on the ami = 0.005 and 
ami =0.01 ensembles all 4 source time slices were binned 
together, except for the 3x3 matrix correlator in the de- 
termination of the T(35) energy. For the latter no bin- 
ning over source time slices was done; instead the error 
estimates from the fits were corrected by 20% upwards 
to be safe. For the ami — 0.02 and ami — 0.03 ensem- 
bles, which have about four times fewer configurations, 
the T(35) state was not computed. There, for the 2x2 
matrix correlators no binning over source time slices was 
performed, again increasing the error estimates by 20% 
upwards instead. For the D-wave correlators, only the 
ID smearing function was included in the fits. Thus, the 
data covariance matrix was small and binning over all 4 
source time slices was used. 



III. RESULTS 

A. Tuning the bare b quark mass 

The bare b quark mass, which is a free parameter in 
the NRQCD action, was tuned non-perturbatively. ft 
was adjusted such that the kinetic mass of the 776(15) 
meson as calculated on the lattice matches the experi- 
mental value of 9.389(5) GeV [n|. The tuning was done 
on the most chiral (ami = 0.005) ensemble of gauge con- 
figurations. 

The kinetic mass was computed from (|11[) . where the 
smallest possible lattice momentum a\p\ = 1 • 2tt/L was 
used. As shown in the next section, the kinetic mass 
is very stable and shows no significant dependence on 
p even for much larger momenta. In order to increase 
statistics the results were averaged over the different pos- 
sibilities for the direction of p. 

The comparison with experiment of course requires 
the knowledge of the lattice spacing, which was de- 
termined as the ratio of the experimentally measured 
T(25) - T(15) mass splitting, 0.56296(40) GeV 0, to 
the dimensionless lattice result. This will be discussed in 
more detail in Sec. IIII CI 

The lattice results for aM kin and the T(25) - T(15) 
splitting at the three different bare quark masses omi = 
2.30, 2.45 and 2.60 are shown in Table [EVl As can be 
seen, the T(25) — T(15) splitting is very insensitive to 
the value of the b quark mass. It is also expected to 
have much smaller lattice discretization errors than the 
IP — 15 splitting as discussed in the next section. 

It was found that in the range considered here, the 



ami, 


aAfkin(?7f)) 


T(2S) - T(1S) 
splitting 


2.30 


4.988(12) 


0.3258(47) 


2.45 


5.281(13) 


0.3242(46) 


2.60 


5.575(13) 


0.3231(54) 



TABLE IV: Results for the tuning of the bare b quark mass 
in lattice units. Errors are statistical/fitting only. 
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FIG. 2: The kinetic mass of the rjt,(lS) meson plotted against 
the bare heavy quark mass. Errors are statistical/fitting only. 
The line shows the average over the bootstrap ensemble of 
linear fit results. 



dependence of the kinetic mass on the bare heavy quark 
mass is described very well by the linear relation 



aMkin = A + B ■ ami,. 



(12) 



A plot of aATkjn as a function of ami is shown in Fig. [2] 
Fits of Eq. (fT2")) with the parameters A and B were per- 
formed on 500 bootstrap samples for the kinetic masses 
at am^ = 2.30, 2.45 and 2.60. The resulting average fit 
parameters were 



A 
B 



0.489(25), 
1.956(11). 



(13) 



To obtain a first result for the lattice spacing of the 
ami — 0.005 ensemble, the T(25) — T(15) mass split- 
ting at ami, = 2.45 was used, giving a" 1 = 1.736(25) 
GeV (the error is statistical/fitting only). Of course the 
6 quark mass was not yet tuned, but given the relative 
independence of the T(25) — T(15) splitting on mj, the 
value of ami, = 2.45 was sufficiently close to the phys- 
ical value. The final results for the lattice spacing ob- 
tained with the correct b quark mass will be presented in 

Sec. Unci 

Using the preliminary result for a -1 , it follows that 
the 77t (15) mass in lattice units must be tuned to be 
a A/km = 5.407(77). Inserting this into (fT2")) and solving 
for a?«6 gives 



am b = 2.514(36). 



(14) 
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n 2 


aM kin (r] b ) 


aC 


c 2 


1 


5.450(17) 


2.5913(84) 


- 


2 


5.450(17) 


2.5912(85) 


1.00003(85) 


3 


5.450(18) 


2.5911(92) 


1.0001(16) 


4 


5.461(22) 


2.597(11) 


0.9981(21) 


5 


5.457(20) 


2.595(10) 


0.9987(24) 


6 


5.452(20) 


2.592(10) 


0.9997(27) 


8 


5.454(22) 


2.593(11) 


0.9993(35) 


9 


5.447(20) 


2.590(10) 


1.0005(35) 


12 


5.445(21) 


2.589(11) 


1.0009(42) 



TABLE V: Kinetic mass, NRQCD energy shift and the square 
of the "speed of light" for various lattice momenta ap = n ■ 
2n/L, calculated on the ami = 0.005 ensemble with amb = 
2.536. 



The error quoted here is statistical/fitting only and is 
dominated by the uncertainty in the lattice result for the 
T(25) - T(15) splitting. 

All remaining calculations were actually performed 
with amb — 2.536. This was an earlier result and the 
fits have been improved slightly since then. However it 
is still inside the range of the new value (fT¥j) . 

For ami, — 2.536 the results were aM^ in — 5.449(13) 
and a" 1 = 1.740(25) GeV. This gives M kin = 9.48(14) 
GeV which is compatible with the experimental result 
of 9.389(5) GeV, confirming the successful tuning of the 
heavy quark mass. 

Note that the lattice NRQCD action can be used 
for both heavy-heavy mesons and heavy-light hadrons. 
Thus, the result for the bare b quark mass obtained here 
will be useful also in future calculations for heavy-light 
hadrons. 



B. Speed of light 

In order to examine how well the lattice data approxi- 
mates the relativistic continuum dispersion relation (fTT)|) , 
the kinetic mass of the ?7t (IS 1 ) meson, defined by (fTT|) . was 
also computed for larger lattice momenta ap = n ■ 2tt/L 
up to n 2 = 12. For these calculations, the local smear- 
ing function T(r) = S r ^o was used at source and sink so 
that multiple lattice momenta can be obtained with lit- 
tle computational cost. For each value of n 2 , the results 
were averaged over the possible directions of the vector 
n, and all components of n were chosen to be less than 
or equal to 2. 

The results are given in Table [V] where also the 
NRQCD energy shift, calculated as 



C 



M kin (p) - E(0) 



and for n 2 > 1 the square of the "speed of light'' 



[E(p) - E(0) + M kinil ] 2 



M 2 



(15) 



(16) 



are shown. In Eq. ([16)) , Mki n ,i denotes the kinetic mass 
calculated with n 2 = 1. In the units used here, one 
should have c 2 = 1. Deviations of c 2 from 1 can be 
caused by discretization errors in the NRQCD, gluon and 
sea quark actions and also by missing higher order rela- 
tivistic corrections in the NRQCD action. The NRQCD 
action is highly improved at tree level, and so the most 
significant errors one expects here are those caused by 
missing radiative corrections. 

As can be seen in the table, in the momentum range 
considered here the kinetic mass shows no significant de- 
pendence on p within the small statistical/fitting errors. 
Correspondingly, c 2 remains compatible with 1, with sta- 
tistical/fitting errors less than 0.5%, indicating that the 
effect of the errors mentioned above is small. 

Analogous calculations for the Y(15) meson have been 
performed in [3] with the same NRQCD action but with 
the Luscher-Weisz gluon and the AsqTad sea quark ac- 
tion. There, the deviation of c 2 from 1 in the same mo- 
mentum range was also found to be compatible with 1 
within statistical errors of less than 1%. 



C. Radial/orbital splittings and the lattice spacing 

The lattice results for the various radial and orbital 
energy splittings are listed in Table IVI1 

Systematic errors are known to be smallest for the spin- 
averaged masses, defined as 



(M) = 



Zj(2J+1)Mj 
Ej(2^+1) 



(17) 



However, in most cases not all of the states entering 
Eq. (|17p are known from experiment. For the 15, 2S 
and 35 masses in this section the J = 1 states (T) are 
considered instead of the spin-averages. Note that the 
J = 5-wave states (rjb) enter the spin-averaged masses 
only with a weight of 1/4, and so the influence of system- 
atic errors in the hyperfine splittings is negligible here. 
For the IP and 2P masses, the spin-averages over the \b 
triplet ( J = 0, 1, 2) states are used. The only experimen- 
tally known Z?-wave state [2l| is T2(l-D) with J = 2, and 
this state is therefore considered here. 

In terms of the NRQCD power counting [l2| . radial 
and orbital energy splittings are of order 0(v 2 ), where 
v is the internal speed of the b quarks inside the heavy- 
heavy meson. For bottomonium one has v 2 w 0.1. The 
NRQCD action in use includes all relativistic corrections 
of order 0(v 4 ) (at tree- level), and hence the missing rela- 
tivistic corrections are of order 0(v 6 ). Naively this leads 
to relativistic errors for the radial and orbital splittings 
of 0(v A ) = 1%. However, as discussed in [14| . for energy 
splittings one has to consider the difference between the 
expectation values of the missing operators for the two 
states. This leads to a reduction of the relativistic errors 
for the 25 - 15 splitting to about 0.5%. 

Additional systematic errors for the NRQCD action 
are due to discretization errors and missing radiative cor- 
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am, = 0.005 


ami = 0.01 


ami = 0.02 


ami = 0.03 


T(2S) - T(1S) 


0.3236(46) 


0.3270(73) 


0.330(18) 


0.327(23) 


T(3S) -T(1S) 


0.517(21) 


0.537(23) 


- 


- 


(Xb(lP))-?(lS) 


0.2589(30) 


0.2572(22) 


0.2628(57) 


0.2613(61) 


( Xb (2P))-T(lS) 


0.478(30) 


0.502(26) 


0.511(39) 


0.516(37) 


( X 6(2P)> - (X6(1P)) 


0.219(29) 


0.245(24) 


0.248(35) 


0.255(33) 


T 2 (li?)-T(15) 


0.4080(46) 


0.4194(42) 


0.417(12) 


0.426(12) 



TABLE VI: Results for the radial and orbital energy splittings in lattice units. Errors are statistical/fitting only. 



ami = 0.005 ami = 0.01 ami = 0.02 ami = 0.03 



a^s-is (GeV) 


1.740(25) (19) 


1.722(38)(19) 


1.708(92)(19) 


1.72(12)(2) 


O-ip-is (GeV) 


1.698(19) (65) 


1.709(15) (65) 


1.673(36) (64) 


1.682(40)(64) 



TABLE VII: Results for the inverse lattice spacing obtained from the T(2S) — T(IS') and (xb(lP)) — T(1S) splittings. The first 
error given is statistical/fitting and the second is an estimate of the systematic errors (relativistic, radiative and discretization) 
due to the NRQCD action. Systematic errors due to the gluon and sea quark actions are not included. 





25-15 


IP - 15 


relativistic 


0.5% 


1.0% 


radiative 


0.5% 


1.7% 


discretization 


0.8% 


3.2% 


total 


1.1% 


3.8% 



TABLE VIII: Estimates of the systematic errors due to the 
lattice NRQCD action for the 2S - IS" radial and IP - IS 
orbital splittings [l4] ]. 



rections (beyond tadpole improvement). Estimates of 
these errors for the 25 — IS* and IP — IS* splittings are 
given in Table IVIIII They are taken to be equal to the 
estimates obtained in [14J for exactly the same lattice 
NRQCD action on the "coarse" MILC gauge configura- 
tions, which have a lattice spacing (a" 1 « 1.6 GeV) very 
similar to the ensembles considered here. The reader is 
referred to [HI and [22j for the details. As can be seen 
in the table, systematic errors are much smaller for the 
25 — 15* splitting compared to the IP — 15 splitting. This 
is due to the smaller difference in the wave functions for 
the 25 and 15 states. The 25 — 15 splitting thus allows 
a more reliable determination of the lattice spacing. 

Note that there are also discretization errors due to 
the gluon and sea quark actions. These are difficult to 
quantify at this stage as only data from one lattice spac- 
ing is available. Gauge configurations with a smaller lat- 
tice spacing are currently being generated by the RBC 
and UKQCD collaborations so that a more systematic 
analysis will become possible in the future. In Q, a 
preliminary error estimate of (oAqcd) 2 ~ 4% for the 
calculations of light hadron properties on the current 
ensembles was given. The calculations performed here 
are different in that the domain wall action only enters 
via the sea quarks. The Iwasaki gluon action [2^, [24| 



is renormalization-group-improved and is therefore ex- 
pected to have a better scaling behavior than the unim- 
proved Wilson action. However, this depends on the ob- 
servable considered; see e.g. [25[ for a scaling study of 
the critical temperature and glueball masses. The sta- 
bility of the "speed of light" demonstrated in Sec. IIIIBI 
provides some evidence for the smallness of the effect of 
gluon discretization errors for bottomonium. 

For reference, the discretization errors in the 25—15 
and IP— 15 splittings on the coarse MILC lattices due to 
the Luscher-Weisz gluon action were estimated in 14] to 
be 0.5% and 1.7%, respectively. These errors are propor- 
tional to the difference in the square of the wave function 
at the origin, which is smaller between the 25 and 15 
states. 

The results for the inverse lattice spacings of the 
four ensembles from both the T(25) — T(15) and the 
(X&(LP)) - T(15) splittings are listed in Table EED For 
the most chiral ensemble the 25 — 15 splitting gives 
a -1 = 1.740(25) s tat(19)syst GeV. No significant depen- 
dence on the sea quark mass can be seen within the 
statistical errors, and therefore no extrapolation was at- 
tempted. For comparison, the RBC and UKQCD collab- 
orations have obtained a -1 = 1.729(28) s t a t in the chiral 
limit, using the f2~ baryon mass f||. This is consistent 
with the results obtained here. 

Next, the lattice spacing determinations from the 25 — 
15 splitting were used to convert the other radial and 
orbital splittings from Table I VII to physical units. The 
results are plotted in Fig. [3l Note that the individual 
results for the lattice spacings of the different ensembles 
were used. 

Overall, good agreement with experiment is seen. The 
dependence on the light sea quark mass is found to be 
weak. This is expected since the typical gluon momenta 
inside the bottomonium are much larger than all the 
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FIG. 3: Radial and orbital energy splittings compared to the 
experimental results (indicated by lines). Errors are statisti- 
cal/fitting only and include the uncertainty in the determina- 
tion of the lattice spacing. The 15 and 25 masses, for which 
no error bars are shown, are not predictions of the lattice cal- 
culation as these states are used to determine the lattice scale 
and the overall energy shift. 



FIG. 4: S-wave hyperfine splittings (energies relative to the 
T(15) and T(25) states, respectively) compared to experi- 
ment. Errors are statistical/fitting only and include the un- 
certainty in the determination of the lattice spacing, which 
enters with a factor of 2. Large systematic errors are expected 
as discussed in the text. 



values for the light quark masses used here. However, 
note that large deviations between lattice results and ex- 
periment were previously seen in quenched simulations 
(rif — 0), so the inclusion of 2+1 flavors of dynamical 
light quarks is in fact very important. A comparison be- 
tween quenched and unquenched results can be found in 



D. Spin-dependent energy splittings 

The spin-dependent energy splittings in bottomonium, 
i.e. the fine and hyperfine structure, are of order 0{v A ) 
and hence any sub-leading corrections are missing in the 
NRQCD action used here. Therefore, the relativistic 
errors in these splittings are expected to be of order 
0(v 2 ) w 10%. The spin-dependent energy splittings also 
receive radiative corrections of order 0(a s ), the strong 
coupling constant at the scale set by the lattice spacing. 
This leads to further systematic errors of the order of 
20%, although tadpole improvement reduces the prob- 
lem. Finally, discretization errors are also expected to be 
larger than for the radial and orbital splittings, especially 
for the S'-wave hyperfine splitting as discussed below. 

The results for the spin-dependent energy splittings 
in lattice units are summarized in Table ILX( where the 
errors given are statistical/fitting only. 



1. S-wave hyperfine structure 

Figure [4] shows a plot of the T(15) — r)b(lS) and 
T(25) —rjb(2S) energy splittings, where the previous lat- 
tice spacing determinations from the 2S — IS splittings 



were used to convert to physical units. 

The errors shown are statistical/fitting only but in- 
clude the uncertainty in the determination of the lattice 
spacing. The latter in fact enters with a factor of 2 here, 
as discussed in [22j, due to the resulting uncertainty in 
the physical heavy quark mass (the hyperfine splitting is 
approximately proportional to the inverse of that mass). 
The statistical error in the IS* hyperfine splitting is then 
dominated by far by this uncertainty, while the 2S hyper- 
fine splitting has an intrinsically higher statistical error 
as the state is radially excited. 

The T(15) — r]b(lS) splitting has recently been mea- 
sured by the BABAR collaboration [n|, who found 
71.4i|?(stat) ± 2.7(syst) MeV. This value is indicated 
in Fig. |U The lattice result in physical units for the 
ami = 0.005 ensemble is 52.5 ± 1.5(stat) MeV, which is 
too small by about 25%, in line with the large systematic 
errors expected. Similarly to the radial and orbital split- 
tings, little dependence on the light sea quark mass is 
seen, which is expected for the same reason as discussed 
there. 

Note that in [T3| and [12] a significant dependence on 
the lattice spacing was found, with the result increasing 
toward finer lattices, indicating that a substantial part of 
the deviation is due to discretization errors. The hyper- 
fine splitting is indeed expected to be sensitive to very 
short distances, as the spin-spin interaction potentials in 
simple models contain a delta function at the origin (see 
e.g. H3). 

Finally, note that in where a relativistic heavy- 
quark action was used, the T(1S) — r)b(lS) splitting on 
the same RBC/UKQCD gauge configurations was found 
to be only 23.7 ± 3.7(stat) MeV, a much larger deviation 
to experiment than found here. 
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am; = 0.005 


ami = 0.01 


ami = 0.02 


ami ~ 0.03 


T(LS) - 




0.03017(14) 


0.03033(16) 


0.03102(36) 


0.03145(38) 


T(2S) - 


776 (25) 


0.0137(30) 


0.0120(48) 


0.013(12) 


0.018(16) 


Xw(lP) 


- <x 6 (i^)> 


-0.0207(20) 


-0.0206(18) 


-0.0231(36) 


-0.0175(70) 




- (Xb(lP)) 


-0.0049(14) 


-0.0027(19) 


-0.0059(22) 


-0.0049(41) 




- (X6(1P)) 


0.0071(11) 


0.0058(12) 


0.0082(17) 


0.0064(29) 


h b (lP) - 


- (xt(iP)) 


-0.0026(18) 


-0.0002(21) 


-0.0014(27) 


-0.0058(42) 


Xbi(lP) 


-Xw(l-P) 


0.0158(18) 


0.0176(25) 


0.0173(40) 


0.0126(77) 


Xw(lP) 




0.0120(23) 


0.0088(31) 


0.0137(38) 


0.0113(68) 


h b (lP) - 


-Xw(lP) 


0.0023(16) 


0.0027(16) 


0.0044(35) 


-0.0009(61) 


Ta(lD) 


-^(1£>) 


0.0011(21) 


-0.0012(18) 


-0.0086(70) 


-0.0050(61) 



TABLE IX: Spin-dependent energy splittings in lattice units. Errors are statistical/fitting only. Large systematic errors are 
expected as discussed in the text. 
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FIG. 5: P-wave spin splittings (energies relative to the spin- 
average of the Xb{lP) states) compared to experiment. Errors 
are statistical/fitting only and include the uncertainty in the 
determination of the lattice spacing. Large systematic errors 
are expected as discussed in the text. 



2. P-wave spin-dependent splittings 

A plot of the IP spin-dependent splittings, converted 
to physical units using the previous 2S — IS lattice spac- 
ing results, is given in Fig. [5l 

It shows the energy differences of the Xboi^P), Xbi(^P), 
Xb2^-P) and hb(lP) states to the spin-average of the 
triplet (xft(lP)). The experimental results [2Cj] for the 
triplet states are also indicated in the plot; the hb states 
have not yet been observed. 

The lattice results are found to be in relatively good 
agreement with experiment, even within the purely sta- 
tistical/fitting errors shown in the plot (those include the 
uncertainty in the lattice spacing). This indicates that 
discretization errors may be smaller than for the .S'-wave 
hyperfine splittings. Note that in simple potential models 
the wave function at the origin is zero (cf. the smearing 
functions in Table |TT| and hence the P-wave spin split- 



tings are expected to be not as sensitive to very short 
distances as the S'-wave hyperfine splittings. 

The result for the experimentally unknown hb(lP) — 
(Xb(lP)) splitting on the ami = 0.005 ensemble is — 4.5 ± 
3.1 MeV, where the error quoted is statistical/fitting only 
and includes the uncertainty from the determination of 
the lattice spacing. 



3. D-wave spin- dependent splittings 

Here, only the T 2 (l-D) — T)b(lD) splitting was calcu- 
lated using the E~~ and representations, as these 
two states do not mix and can be computed from the 
same heavy-quark propagators. 

The lattice results for the different ensembles are listed 
in Table llXl On the am; = 0.005 ensemble, the splitting 
in physical units is found to be 1.8 ± 3.7 MeV where 
the error given is statistical/fitting only and includes the 
uncertainty from the determination of the lattice spacing. 
No experimental results are available. 



IV. CONCLUSION 

In this paper, a comprehensive calculation of the bot- 
tomonium spectrum with improved lattice NRQCD on 
the RBC/UKQCD 24 3 x 64, L s = 16 gauge configura- 
tions with 2+1 flavors of dynamical domain wall fcrmions 
was presented. The results are similar to those ob- 
tained in with the same heavy-quark action from the 
"coarse" MILC gauge configurations, which use the Asq- 
Tad fermion action and the Liischer-Weisz gluon action. 
In particular, good agreement with experiment was found 
for the radial and orbital energy splittings, for which sys- 
tematic errors due to the NRQCD action are small. Fur- 
thermore, no significant deviations of the "speed of light" 
from 1 in the rjb(lS) dispersion relation were found within 
the small statistical errors. The calculations in this work 
provide further evidence for the good properties of the 
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domain wall and Iwasaki actions employed by the RBC 
and UKQCD collaborations. By comparing the 2S — IS 
radial energy splitting to experiment, independent deter- 
minations of the lattice spacings were performed, giving 
a -1 = 1.740(25) sta t(19) S y S t GeV for the most chiral en- 
semble. 

The results for the fine and hyperfine structure are 
expected to have larger systematic errors due to miss- 
ing radiative and relativistic corrections in the NRQCD 
action as well as discretization errors. Nevertheless, rela- 
tively good agreement with experiment was seen for the 
P-wave fine structure, and the deviation to experiment in 
the 15 hyperfine splitting was found to be much smaller 
than for the previous result obtained from a relativistic 
heavy quark action in jioj ]. 

The calculations presented here are only for one lat- 
tice spacing. A more systematic analysis of discretiza- 
tion errors will be performed once new ensembles with a 
finer lattice spacing are made available by the RBC and 
UKQCD collaborations. 

Having obtained the bottomonium spectrum and the 
bare heavy quark mass in this work, the next step is 
to perform calculations for heavy-light systems. Results 
for the spectrum of heavy-light baryons and mesons with 
domain wall valence quarks but with static heavy quarks 
were recently presented in [28]. Heavy- light computa- 



tions with NRQCD heavy- and domain wall light valence 
quarks on the same RBC/UKQCD gauge field ensembles 
are currently underway. 
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